library(sf)
library(stars)
library(caret)
library(CAST)
library(sen2r)
library(tmap)
library(knitr)
We will first fix the random number seed, to get identical results for procedures that involve randomness. Remove this command if you want the random effect in outcomes.
set.seed(131)
trainsites <- st_read("../data/trainingsites_muenster.gpkg")
## Reading layer `trainingsites_muenster' from data source `/home/hanna/Documents/Github/HannaMeyer/IGARSS_21/data/trainingsites_muenster.gpkg' using driver `GPKG'
## Simple feature collection with 29 features and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 400682 ymin: 5754417 xmax: 408327.9 ymax: 5760113
## projected CRS: WGS 84 / UTM zone 32N
list_safe <- s2_list(
spatial_extent = st_as_sfc(st_bbox(trainsites)),
tile = "32ULC",
level = "L1C",
time_interval = as.Date(c("2020-04-26", "2020-04-28"))) # ursprüngliche szene: c("2019-04-17", "2019-04-19")
#s2_download(list_safe[1], outdir="../data/")
sen_id <- names(list_safe[1])
bands <- c("B04", "B03", "B02", "B08")#, "B05", "B06", "B07", "B8A", "B11", "B12")
s2 <- list.files(paste0("../data/",sen_id),
recursive=TRUE,
full.names = TRUE,
pattern=".jp2")
# match band name to file name:
m <- match(paste0(bands,".jp2"), substr(s2,nchar(s2)-6,nchar(s2)))
s2 <- s2[m]
sen <- read_stars(s2, proxy = TRUE, NA_value = 0) %>%
setNames(bands)
pts <- st_as_sf(st_sample(trainsites, 200, "random"))
pts <- st_intersection(pts,st_make_valid(trainsites))
trainDat <- st_extract(sen, pts) %>%
st_intersection(pts)%>%
data.frame()
trainDat <- data.frame(trainDat,pts)
trainDat$Label <- as.factor(pts$Label)
ind <- CreateSpacetimeFolds(trainDat,spacevar="PolygonID",class="Label",k=3)
ctrl <- trainControl(method="cv",index=ind$index,indexOut = ind$indexOut,savePredictions = TRUE)
model <- train(trainDat[,attributes(sen)$names],
trainDat$Label,
tuneLength = length(attributes(sen)$names)-1,
trControl = ctrl,
method="rf",ntree=200)
prediction <- predict(sen, model)
confusionMatrix(model$pred$pred[model$pred$mtry==model$bestTune$mtry],
model$pred$obs[model$pred$mtry==model$bestTune$mtry])
## Confusion Matrix and Statistics
##
## Reference
## Prediction Fallow field Grassland Inland water Mixed forest Planted field
## Fallow field 0 0 0 1 19
## Grassland 0 8 0 0 2
## Inland water 0 0 18 0 0
## Mixed forest 1 0 0 46 1
## Planted field 25 1 0 0 3
## Urban 2 0 3 0 13
## Reference
## Prediction Urban
## Fallow field 0
## Grassland 0
## Inland water 0
## Mixed forest 0
## Planted field 3
## Urban 54
##
## Overall Statistics
##
## Accuracy : 0.645
## 95% CI : (0.5744, 0.7112)
## No Information Rate : 0.285
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.5477
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: Fallow field Class: Grassland Class: Inland water
## Sensitivity 0.0000 0.8889 0.8571
## Specificity 0.8837 0.9895 1.0000
## Pos Pred Value 0.0000 0.8000 1.0000
## Neg Pred Value 0.8444 0.9947 0.9835
## Prevalence 0.1400 0.0450 0.1050
## Detection Rate 0.0000 0.0400 0.0900
## Detection Prevalence 0.1000 0.0500 0.0900
## Balanced Accuracy 0.4419 0.9392 0.9286
## Class: Mixed forest Class: Planted field Class: Urban
## Sensitivity 0.9787 0.07895 0.9474
## Specificity 0.9869 0.82099 0.8741
## Pos Pred Value 0.9583 0.09375 0.7500
## Neg Pred Value 0.9934 0.79167 0.9766
## Prevalence 0.2350 0.19000 0.2850
## Detection Rate 0.2300 0.01500 0.2700
## Detection Prevalence 0.2400 0.16000 0.3600
## Balanced Accuracy 0.9828 0.44997 0.9107
AOA <- merge(sen)%>%
split() %>%
st_as_stars(downsample=c(10,10))%>%
setNames(bands)%>%
aoa(model)
ext_muenster <- st_buffer( st_as_sfc(st_bbox(trainsites)),1500)
png("../figures/rgb_all.png", width=15, height=15,units = "cm",res=300)
plot(st_rgb(merge(st_as_stars(sen,downsample=c(10,10)))[,,,c(1,2,3)], stretch = TRUE, probs = c(.01, .99)),main=NULL,reset=FALSE)
plot(ext_muenster,add=TRUE,lwd=3,border="red")
invisible(dev.off())
cols <- rev(c("red","lightgreen","forestgreen","blue","green","beige"))
prediction_st <- st_as_stars(prediction, downsample = c(10,10))
map1 <- tm_shape(prediction_st, raster.downsample = FALSE) +
tm_raster(palette = cols,title = "LUC")+
tm_scale_bar(bg.color="white")+
tm_layout(legend.position = c("left","top"),
legend.bg.color = "white",
legend.bg.alpha = 0.6,
legend.text.size = 1.5,
legend.title.size = 1.5)
prediction_st[AOA[2] == 0] <- NA
map1_AOA <- tm_shape(prediction_st, raster.downsample = FALSE) +
tm_raster(palette = cols,title = "LUC")+
tm_scale_bar(bg.color="white")+
tm_layout(legend.position = c("left","top"),
legend.bg.color = "white",
legend.bg.alpha = 0.6,
bg.color="black",
legend.text.size = 1.5,
legend.title.size = 1.5)+
tm_add_legend(type = "fill",
col="black",
labels = "Outside AOA")
tmap_save(map1, "../figures/LUC_map.png")
tmap_save(map1_AOA, "../figures/LUC_map_AOA.png")
sen_crop <- st_crop(sen,st_bbox(ext_muenster))%>%
st_as_stars()%>%
merge()
prediction_crop <- st_crop(prediction,st_bbox(ext_muenster))
AOA_ms <- aoa(sen_crop, model)
prediction_st_crop <- st_as_stars(prediction_crop)
png("../figures/rgb_ms.png", width=15, height=11,units = "cm",res=300)
plot(st_rgb(sen_crop[,,,c(1,2,3)], stretch = TRUE, probs = c(.01, .99)),main=NULL,reset = FALSE)
plot((st_geometry(trainsites)),add=TRUE,border="red",lwd=3)
invisible(dev.off())
map1_ms <-tm_shape(prediction_st_crop, raster.downsample = FALSE) +
tm_raster(palette = cols,title = "LUC")+
tm_scale_bar(bg.color="white")+
tm_layout(legend.position = c("left","top"),
legend.bg.color = "white",
legend.bg.alpha = 0.6,
legend.text.size = 1.5,
legend.title.size = 1.5)
prediction_st_crop[AOA_ms[2] == 0] <- NA
map1_AOA_ms <-tm_shape(prediction_st_crop, raster.downsample = FALSE) +
tm_raster(palette = cols,title = "LUC")+
tm_scale_bar(bg.color="white")+
tm_layout(legend.position = c("left","top"),
legend.bg.color = "white",
legend.bg.alpha = 0.6,
bg.color="black",
legend.text.size = 1.5,
legend.title.size = 1.5)+
tm_add_legend(type = "fill",
col="black",
labels = "Outside AOA")
tmap_save(map1_ms, "../figures/LUC_map_ms.png")
tmap_save(map1_AOA_ms, "../figures/LUC_map_ms_AOA.png")
Run this code only if you want to write the results. It takes a while!
#write_stars(AOA,"../data/AOA.tif",layer=2)
#write_stars(prediction,"../data/prediction.tif",
# chunk_size=c(dim(prediction)[1],
# floor(2.5e+06/dim(prediction)[1])))